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Abstract 

We perform accurate investigation of stability of localized vortices in an effectively two- 
dimensional ( "pancake-shaped" ) trapped BEC with negative scattering length. The analysis com- 
bines computation of the stability eigenvalues and direct simulations. The states with vorticity 

( 9=1 ^ 

iS* = 1 are stable in a third of their existence region, < < (l/3)iVAiax , where is the number of 
(5=1) 

atoms, and A^max is the corresponding collapse threshold. Stable vortices easily self-trap from ar- 

bitrary initial configurations with embedded vorticity. In an adjacent interval, (l/3)A^max < N < 
(5=1) 

0.43A^max , the unstable vortex periodically splits in two fragments and recombines. At > 

(5=1) 

0.43A^max , the fragments do not recombine, as each one collapses by itself. The results are 
compared with those in the full 3D Gross-Pitaevskii equation. In a moderately anisotropic 3D con- 
figuration, with the aspect ratio Vw, the stability interval of the S = 1 vortices occupies ^ 40% 
of their existence region, hence the 2D limit provides for a reasonable approximation in this case. 
For the isotropic 3D configuration, the stability interval expands to 65% of the existence domain. 
Overall, the vorticity heightens the actual collapse threshold by a factor of up to 2. All vortices 
with S >2 are unstable. 

PACS numbers: 03.75.Lm,03.65.Ge,05.45.Yv 
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I. INTRODUCTION 

Vortices are fundamental dynamical objects in Bose- Einstein condensates (BECs) 
They can be easily created in large numbers, forming vortex lattices Q]. A review of basic 
results for the BEC vortices can be found in Ref. In most works, vortices were studied 
in self-repulsive BECs, with a positive scattering length characterizing collisions between 
atoms. In the experiment, vortices are created stirring the condensate by a properly designed 
laser beam ^, or by imprintingan appropriate phase pattern onto a condensate trapped in a 
Joffe-Pritchard magnetic trap . Recently, more complex vortical structures were predicted 
in repulsive condensates, such as vortex dipoles and quadrupoles , and patterns in the 
form of globally linked topological defects 

Vortices in repulsive BECs may be regarded as dark solitons on a finite background, 
similar to optical vortices in a self-defocusing dielectric medium j^. A challenging issue 
is the study of vortex solitons (completely localized objects with embedded vorticity) in 
attractive BECs, with a negative scattering length j^. In particular, attractive interaction 



occurs in the ''Li condensate where quasi-one-dimensional solitons were created 

The most natural setting for vortex solitons is a "pancake" condensate, strongly confined 
in one direction (z) and weakly trapped in the radial direction (r) in the transverse plane. 
In the experiment, the "pancake" is created by a superposition of a tight optical trap, with 
a confinement thickness a^, and a loose radial magnetic trap with a trapping frequency 
ilr in the transverse plane jl^. The theoretical description of the "pancake" configuration 
assumes that the underlying three-dimensional (3D) Gross-Pitaevskii equation (GPE) may 
be reduced to its 2D counterpart. The condition necessary for the reduction is that the 
squared harmonic-oscillator length determined by the radial trapping frequency Qr, 

{a^^^'^Y = h/ {mQr) (1) 

(m is the atomic mass), must be much larger than a^, i.e., 

-C Tc^h/ (mal) . (2) 

For Li atoms trapped in the diffraction-limited gap between two parallel strongly repulsive 
(blue-detuned) light sheets, with — 2 fim, this means Qr ^ 10 KHz, which is readily 
met in the experiment, where confining frequencies are measured in tens of Hz. Condition 
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(j21) also rule out bending instability of vortex cores in the pancake geometry (in the full 
3D case, the instability results in decay of straight vortices into vortex rings, in repulsive 
condensates Q]), as is much smaller than the vertical length necessary for the bending. 

The condition (0) for the applicability of the 2D approximation can be cast in another 
form, which sets a limit on the number of atoms M in the pancake-shaped condensate. As 
shown below, this form is 

u < j: , , (3) 

2 /im as above, 
10 /zm 



h 

mVLr\a\az ' 

where a is the (negative) scattering length. Taking the same value 
together with experimentally relevant Vt^ = 10 Hz [in this case, Eq. yields a^^°^ 



-0.1 nm, which is typical for soliton experiments in ''Li 
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ll|,Eq. 



for Li atoms], and a 
© yields ^^ < 10^. 

Generally, vortex solitons are unstable against collapse (catastrophic self- compression of 
the wave function), and also against azimuthal perturbations which break the axial symmetry 
of the solitons. In fact, condition (jS)) provides for the stability of the quasi-2D solitons 
against the 3D collapse. Within the framework of the 2D description, both the ordinary 
(zero-vorticity) solitons and their vortical counterparts can be stabilized by a square-shaped 
optical lattice (OL) 



U, 3| • The same OL can also support weakly localized vortex solitons 
of the gap type in a repulsive BEG Q, [l^; a quasi- ID lattice stabilizes ordinary solitons, 



but not vortices, in the 2D model with self-attraction 
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3- 



In attractive condensates. 



ordinary 3D solitons can be stabilized by low- dimensional OLs [viz., a quasi-2D lattice 
IgI 18 . Iigj l. or an axially symmetric one 2Q|], as well as by the full 3D lattice 21 1 (for a 
recent review of multidimensional solitons in optics and BEG, see Ref. j^j)- In addition, 2D 
zero-vorticity solitons can be stabilized by means of the Feshbach-resonance management 
[2^, which amounts to time-periodic modulation of the scattering length between positive 



and negative values 



24| . However, the latter mechanism cannot stabilize vortex solitons 



26|. 



25| . although it may inhibit the development of their instability 

While vortex solitons of the 2D GPE with a negative scattering length are strongly 
unstable in free space, they can be stabilized by the radial trappi ng p otential, (m/2)fi^r^. 
This issue was considered in an early work [2^ and more recently jsJSj- ^^^^ known 
that the radial trap stabilizes zero-vorticity {S = 0) states against collapse, provided that 
the norm of the solution (which is proportional to the number of atoms in the condensate) 
is limited from above. In the 3D case with the isotropic radial trap, the maximum norm 



depends on the corresponding trapping frequency Qr- This maximum norm was found 
by means of the variational approximation and from numerical computations ■ In 

contrast to that, the maximum 2D norm, A^i^x^'', does not depend on fi^, being exactly 

n 

equal to the norm of the corresponding Townes soliton in the free 2D space 30] . Similarly, 
the maximum value of the norm for the trapped fundamental vortex (one with S = 1), which 
does not lead to collapse, depends on VLr in the 3D case (the shape of the 3D trap which 



maximizes the collapse threshold under additional conditions was studied in Ref. [3l|), but 
not in the 2D situation, where it coincides with the known collapse threshold, A^max , of the 
fundamental vortex soliton in free space j^^]- Formally, the vorticity gives rise to a dramatic 
increase of the 2D collapse threshold: Nmia^'' ~ ANma^^^ ■ However, we will demonstrate that 
the actual increase of the 2D stability threshold, A^thr? in the fundamental vortex is a more 
modest effect, amounting to A^thr — 2A^]^x°''; in the remaining part of the existence interval 
for the vortex, A^thr < N < A^max , the azimuthal instability splits the vortex into a pair of 
fragments, which then collapse intrinsically. 

It is relevant to mention that, in various models of nonlinear optics (see Refs. 133| and brief 
), vortex solitons, while being stable against radial perturbations, are easily split 
into fragments by an azimuthal instability. Nevertheless, azimuthally stable vortical solitons 
have been predicted too, in models with competing 3^ or nonlocal nonlinearities js^, and 
in defocusing media with an imprinted cylindrical lattice 3^ . 

In this work, we aim to investigate the stability of vortices in the 2D GPE equation 
with self-attraction and weak radial trapping potential, through computation of a full set of 
the corresponding stability eigenvalues. The results will be verified by direct simulations. 
We conclude that all vortices with S >2 are unstable (as conjectured in Ref. Q])' while 
the fundamental ones {S = 1) are completely stable for < < A^^cr ~ (l/3)A^max [in 
Ref. Q, it was conjectured that N^^ ^ (l/4)A/'^i;^^]. At A^ = A^'^, the S" = 1 vortex 
is destabilized by symmetry- breaking oscillatory perturbations. Further, in an adjacent 

(S = l) (S=D 

interval, (l/3)A^max < A^ < 0.43A^max , the evolution of the unstable vortex does not 
lead to collapse; instead, it features nearly periodic splittings in two fragments followed by 
their recombinations back into the vortex. This regular dynamical regime was not known 
before. In the remaining part of the existence interval, 0.43A^max < N < A^max , fragments 
produced by the splitting of the vortex do not recombine, but rather blow up intrinsically, 
as their individual norm exceeds the collapse threshold for the Townes soliton. 
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We will also verify that the results obtained within the framework of the 2D model are 
valid indeed in the full 3D GPE with the corresponding anisotropic trap. In particular, in 
a moderately anisotropic 3D model, with the radial-confinement length a^'^'' [see Eq. (P)] 
larger by a factor of VTO than its counterpart, a^^"\ in the orthogonal (tightly confined) 
direction, the relative size of the stability area inside the existence region for 3D solitons 
with the embedded vorticity S" = 1 is ~ 0.4, which is to be compared to the above-mentioned 
relative size of the stability region, 1/3, in the 2D limit. Approaching the isotropic limit, 
a^^"^ /a^"^ — > 1, the stability area expands, attaining the relative size ~ 0.65. 

The paper is organized as follows. In Section 2, we formulate the 2D model, briefly 
recapitulating its derivation from the 3D GPE. In Section 3, we present basic results for the 
2D model: shapes of the vortex states (actually, we find them by continuation of exact states 
with a definite value of the angular momentum in the linear 2D Schrodinger equation), and 
eigenvalues that determine their stability. In Section 4, we verify the predicted stability 
conditions by direct simulations. In Section 5, we consider the full 3D model with the 
anisotropic trap, and Section 5 concludes the paper. 



II. THE MODEL 



We assume a self-attractive condensate which is loaded in a nonrotating loose trap in 
the [x, y) plane, and tightly confined in the z direction. The corresponding 3D GPE for the 
single-atom wave function \1/ is 38]: 

where is the 3D Laplacian and a is the negative s-scattering length accounting for the 
attraction between atoms. In the pancake configuration, the confinement frequency VLz is 
assumed to be much larger than Vtr- If the tight confinement in the z-direction is provided 
by parallel light sheets repelling the atoms, the corresponding potential term, instead of 
mf22^^/2, is a deep rectangular potential well in the region of < ^ < a^. 

The 3D equation can be reduced to an effective GPE in two dimensions, provided that 
the transverse quantum pressure is much stronger than self-attraction. In other words, the 
energy of the atom in the ground state of the vertical trap, if the latter is taken as the deep 
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potential well of width Oz, 

Eo = {nhf / {2mal) , (5) 

must be much larger than the contribution of the attraction to the atomic chemical poten- 
tial, Afi ~ —2TT?i^{a/m)n, where n = |\E'|^ is the 3D atom density [if the tight vertical 
confinement is provided by large Qz in Eq- 0? then the respective harmonic-oscillator 
length, a^^°^ = ^Jfij (mfi^), taken as per Eq. ((H), defines the effective confinement size as 

{ho)i 

For the pancake-shaped configuration in the weak radial trap, the density is determined, 
together with a characteristic radial size R of the pancake itself, by the condition of the 
balance between the radial quantum pressure, radial trapping force, and nonlinear self- 
attraction, which yields 

V mfi, " ^ ' ^Tiah ■ ^ ' 

The substitution of expressions ^ and (jH)) in the above condition ^ I^A^I? which 
implies the tight transverse confinement, leads to Eq. (j21) that determines the range 
of the radial trapping frequencies where the 2D description is applicable. As a conse- 
quence of this condition, the aspect ratio of the pancake configuration is always large, 
TiR/az ~ a'^^Tij {mQr) ^ 1. A limitation on the number of atoms in the pancake state 
can be derived by noting that it is determined by the product of the density, as given by Eq. 
(jni), and the effective volume, ~ vri^^a^, which leads to Eq. (jH)). If A/" exceeds this limit, the 
condensate will start a nearly-2D collapse, which may later go over into a stronger collapse 
in three dimensions. 

The final derivation of the effective 2D equation follows a well-known route: the 3D wave 
function is factorized, 

^ix,y,z,t) = 'ip{x,y,t)(f)oiz), (7) 

where 0o = sinlnz/az) is the ground-state wave function of the transverse state between 
two hard walls separated by the distance a^. After that, aver agin g of the 3D equation (jlj 
in the z direction leads to a normalized two-dimensional GPE 13911 . 
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i/-, (8) 



where g = Sna (or g = 2^/2Tca for the case of the parabolic trap in the vertical direction), 
the transverse-ground-state energy Eq was subtracted from the chemical potential, while Qr 



are 



and t actually stand for mVLr/Ti and {h/m)t. The variables x, y and IV'P in Eq. i 
measured in the same units as in Eq. We then redefine them too, hy ip ^ y^Jg\/Tl^'i/j, 
{x,y) \/ri^{x,y), and additionally rescale time by t ^ Qrt, to finally set = —g = 1 
(unless g = 0). 

Equation (jHl) conserves the 2D norm, = j J \tlj{x,y,t)f dxdy, and energy 



E 




dx 



+ 



dil) 



dy 



-^gli'l' + l^lix' + y'Ml' 



dxdy. 



(9) 



A straightforward consequence of Eq. (jH)) is the relationship E = fiN + TT\g\ rR^dr, 
which is used to cross-check the accuracy of numerical results, see below. Further, the above 
rescalings and factorization (|7j) result in the following relation between the 2D norm and 
the number of atoms. A/" (which is defined as the norm of the 3D wave function in original 
units). 



^^ 



-N = 



(10) 



2\g\ 67r|a| 

where the second equality follows from the above relation g = Stto, if the vertical confinement 
is provided by the repelling light sheets. If, instead, a tight parabolic potential is used to 
confine the condensate in the 2;-direction, the coefficient in front of the last term in Eq. (|iup 
is replaced by ai'^°'*/ (2-\/2|a|), where ai^°'* is the respective linear-oscillator length as defined 
above. 

The final normalized form of the 3D GPE, which is also considered below, is 



+ 



(11) 



2 \dx^ dy"^ dz^ ) ' 2 
Here the nonlinearity coefficient is again scaled to be —1, and = Q'l/Q'^ measures the 
relative strength of the tight confinement in the vertical direction, the nearly-2D case cor- 
responding to ^ 1. The norm of the 3D solutions, N = J J J {"^{xjy, z,t)f dxdydz, is 
related to the number of atoms as follows: 

M = \ — — = ^— AT, 12) 

V milr 87r|a| Svrlal 

where a^°^ is the harmonic-oscillator length in the radial direction, defined as per Eq. 
The conserved energy of the 3D condensate is 



E 







(9^ 
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III 




dx 


+ 


dy 


+ 


dz 



1^1 



dxdydz. 
(13) 



Note that, after Qr and —g were normalized to be 1, the two-dimensional GPE (jSJ features 
no free parameters, while its 3D counterpart, Eq. (|TT|. contains one free coefficient, Q^. 

III. TWO-DIMENSIONAL VORTEX SOLUTIONS AND THEIR LINEAR STA- 
BILITY 

Stationary vortex solutions to Eq. (jH)) are looked for as 

= R{r) exp {iSe - ifit) , (14) 
with an integer vorticity S, real chemical potential /i, and function R{r) obeying the equation 

R" + r~^R! + - - r^) i? - 2gR^ = 0, (15) 

where the prime stands for d/dr. In the linear case, g = 0, Eq. (fT3j) is tantamount to the 
usual stationary quantum-mechanical Schrodinger equation for the 2D harmonic oscillator, 
which gives rise to an infinite set of discrete eigenvalues, fi = j + k + 1, with j,k = 0, 1, 2, ... 
. In the Cartesian coordinates, the corresponding eigenfunctions are (recall we have set 
Q = l) 

^jk{x,y,t) = e-^(^+'=+^)*$,(x)$fc(y), (16) 

where $j(x) and ^k{x) are stationary wave functions of the ID harmonic oscillator associ- 
ated with the energy eigenvalues j + 1/2 and A; + 1/2, respectively. The set of states (fT^ 
is degenerate, as all the states corresponding to a fixed value of j ' + A; appertain to the 
same eigenvalue. Making use of this degeneracy, the eigenfunctions may be rearranged into 
combinations with a definite value of the vorticity, which is 

S = j + k = i2-l. (17) 

In particular, the combination 

4^=') = ^^o{x, y, t) + #01 (a;, y,t)=r exp {-2tt + i6 ~ /2) (18) 

is the wave function with the lowest value of the chemical potential for 5* = 1, ii = j + k+l = 
2, and 

4"^=') = ^2o(a:, y, t) - V'o2(a;, y, t) + 2#ii(a;, y, t) = exp (-3^t + liB - jl) (19) 
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is the wave function with the lowest chemical potential, /i = j + A; + l = 3, for S = 2. 

Nonlinear solutions for vortices can be obtained by a continuation method, starting at 
(yf = with expressions (fTH|) and (fT^ . Note that, pursuant to Eq. (fT^ . the shift of /i due 
to the self-attractive nonlinearity (with 51 < 0) is negative, hence from the comparison with 
the linear result (fTTj) it follows that 

^l{N ^ 0) < ^l{N = Q) = /i^ax = S + l. (20) 

To specify solutions to Eq. (|15p in the nonlinear case (recall we set g = —1), Eq. (|15p 
was solved numerically. The dependences fi{N) obtained from the numerical solutions for 
S* = 0, 1 and 2 are plotted in Fig. 1(a) (in a part, these dependences are tantamount to those 
reported earlier in Ref. j3]). Naturally, fi approaches the linear value //max for N ^ 0, see 
Eq. pO|) . while for large negative fi the norm of the S = state asymptotically approaches 
the value A^i^x^'* ~ 5.85, which is the above-mentioned collapse threshold corresponding 
to the Townes soliton |30]. Similarly, for the states with 5* > 1 the asymptotic values 
N{fi = —00) are the above-mentioned ones, A^max, which determine the (formal) collapse 
threshold for the vortices. 

To study stability of the stationary solutions, we search for a perturbed solution of Eq. 
dHl) as 

ip{x, y, t) = [R{r) + u{r) exp(At + iL9) + v*{r) exp(A*t - iLO)] exp {iS9 - ifxt) , (21) 

where {u, v) and A are eigenmodes and the instability growth rate corresponding to an 
integer azimuthal index L of the perturbation. The linearization around the stationary 
solution leads to equations 

iXu + ^ [u" + r'^u' - (5* + Lfr^'^u] + fiu + R^{v + 2u) - K'^u = 0, 
-iXv + ^ [v" + r~^v' - (S - L)V-\] +fiv + R^{u + 2v) - Vt; = 0, (22) 

which were solved numerically, with boundary conditions demanding that u{r) and v{r) 
decay exponentially at r — 00, and decay as r'"^^^' at r — 0. 

The most important result following from the numerical computation of eigenvalues A is 
that the vortex soliton with S" = 1 is stable, i.e., real parts of all the eigenvalues remain 
equal to zero, only in a finite interval of values of the chemical potential adjacent to the 
linear limit [see Eq. (|^]. fJ.nia.x{S = 1) = 2 > > ~ 1.276, and outside this interval, 
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i.e., for — oo < jj, < ficr, the states with 5* = 1 are unstable. At the critical point (/i = /icr), 
the corresp 
the region 



the corresponding critical norm is A^ir^^ = 7.79. Thus, the vortices with S* = 1 are stable in 



< < N^^^^ = 7.79, (23) 
and in the remaining two-thirds of their existence interval, i.e., 

7.79 < AT < A'£=i) ^ 24.1, (24) 

they are unstable. The conclusion that "less nonlinear" solutions, i.e., ones with smaller A^, 
are more stable is quite natural, as they continue the stable linear solutions which correspond 
toA^ = 4^. 

The energy [see Eq. Q] at = /icr is Ecr = 12.913. The energy as a function of A^ is 
displayed, for the S = 1 soliton family, in Fig. 1(b). It is noteworthy that, while the fi{N) 
dependence in Fig. 1(a) is monotonic, its E{N) counterpart is not. 

To further characterize the (in) stability of the vortices with S = 1 and 2, in Fig. 2(a) we 

(S) 

show the largest real part of the eigenvalues A]^ as a function of jj, for different values of the 
perturbation azimuthal index L, see Eq. ()21|) . Extensive numerical calculations demonstrate 
that the perturbation mode that actually determines the stability border always has L = 2. 
In particular. Re ^A^^^''^ is positive in the entire existence domain of the 5* = 2 vortices, 
up to the linear limit, firna.x{S = 2) = 3 [see Eq. thus making the vortices with S = 2 

completely unstable and confirming a conjecture put forward in Ref. 0|. All vortices with 
S* > 2 are unstable too. 

For the fundamental vortex solitons with 5 = 1, the numerical calculation of the eigen- 

( S=l) 

values was performed for values of the azimuthal index up to L = 6, in a broad range 

of values of the soliton's chemical potential, —20 < fi < 2. It was found that the maximum 
value of Re ^A^^~^^ j is nonzero for L = 1, 2, 3, and this value is zero (up to the numerical 
accuracy) for higher values of the azimuthal index, L = 4,5, and 6. To further check the 
accuracy of the numerical results, we have compared the eigenvalues produced by the nu- 
merical code running on 400 and 800 discretization points, respectively, and concluded that 
four significant digits are identical in both cases. 

The bifurcation responsible for the destabilization of the S = 1 vortex at yU = ficr is 
identified by results displa^d in Fig. 2(b). As seen from the figure, the bifurcation is of the 
Hamiltonian-Hopf type j4l|, i.e., it is accounted for by a quartet of complex eigenvalues, 
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(A, A*, —A, —A*), cf. earlier reported results for solitons in scalar |4^|43| and vectorial [44 . 
models. The unstable quartet is generated by a collision of two pairs of imaginary (stable) 
eigenvalues, which is typical to the Hamiltonian-Hopf bifurcation. 

IV. DIRECT SIMULATIONS OF TWO-DIMENSIONAL STABLE AND UNSTA- 
BLE VORTICES 

To verify the above predictions concerning the stability of the vortex solitons in the 2D 
approximation, we have performed direct simulations of the perturbed evolution of vortex 
solitons in Eq. (jHJ. The simulations were carried out by means of the Crank-Nicholson 
discretization scheme. The system of the corresponding nonlinear finite-difference equations 
was first solved by means of the Picard iteration method 3| , and the resulting linear system 
was then handled using the Gauss-Seidel iterative scheme. For good convergence we typically 
needed five Picard iterations and eight Gauss-Seidel iterations. We employed a grid with 
600 X 600 points (typical stepsizes were Ax = Ay = 0.025 and At = 0.0008). 

To test the stability of the S = 1 vortex, random noise was added to it at the initial 
moment - typically, with a 5% relative amplitude. Figure 3 shows how a vortex with S = 1, 
that was predicted above to be stable against small perturbations, completely recovers after 
the addition of the initial perturbations. 

To further explore the robustness of the vortex which is predicted to be stable, in Fig. 4 
we display its self-trapping from an arbitrarily chosen Gaussian input with a nested vortex, 
Uq = Ar exp (— ar^ + i6), that has A = 1 and a = 1. The shape of the input is far from the 
exact stable vortex soliton. In this case, we used a grid with 391 x 391 points and stepsizes 
As = Ay = 0.028, At = 0.001. Note that, at the values of corresponding to the definitely 
stable solitons displayed in Figs. 3 and 4, a rather inaccurate stability condition proposed 
in Ref. [2^ for the S = 1 vortices would imply instability. 

For the S = 1 vortices which are predicted above to be unstable, we have found that, in 
the interval of 

7.79 < < 10.30, (25) 

which is adjacent to stability region and occupies ~ 1/8 of the entire existence region 
of the vortices [cf. Eq. (j24|) ]. initial perturbations initiate regular quasi-periodic evolution, 
which features recurrent splittings of the soliton into two segments and their recombina- 
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tions. The regular character of the evolution in this case is also manifest in nearly periodic 
oscillations of the soliton's amplitude, see typical examples in Figs. 5 and 6. This dynamical 
regime (which, to the best of our knowledge, has never been reported before) can be easily 
understood, as, unlike the known scenario of the instability development for vortex solitons 
in nonlinear-optical models in free space, where the splinters emerging after the breakup of 
the vortex soliton separate, moving away in tangential directions, in the present setting they 
cannot do it, being confined by the parabolic trap. We stress that Eq. ()25l ) shows that the 
introduction of the vorticity makes it possible to heighten the actual collapse threshold for 
solitons in the trapped self-attractive condensate by a factor of two: from Nma^^^ = 5.85 for 
the spinless solitons to ~ 11 for their vortical counterparts. 

For larger values of the norm, 10.3 < N < A^max ~ 24.1, the evolution of the unstable 
vortex again starts with its splitting into two fragments; however, in this case they do not 
recombine, but rather quickly blow up (not shown here). The fact that the corresponding 
collapse threshold for the S = 1 vortex solitons, found at = 10.3, is much smaller than the 
above-mentioned formal threshold, iVmax ~ 24.1, can be readily explained. Indeed, it has 
been shown above that the azimuthal instability sets in earlier than the radial instability 
that would directly lead to collapse. If the norm of each of the zero-vorticity fragments, 
into which the vortex is broken by the azimuthal instability, is close enough to the collapse 
threshold for the 5 = solitons, i.e., A^i^x"'* = 5.85, the fragments start collapsing and 
therefore fail to recombine into the vortex, unlike the situation with the weak instability in 
interval ((2S1)- This explanation is consistent with the fact that the actual collapse threshold 
for the vortex, = 10.3 , is fairly close to 2Nma^^^ ■ 

As concerns the S* = 2 solitons which were predicted above to be unstable at all finite 
values of A^, the simulations (not shown here) demonstrate their irregular (chaotic) evolution 
for < A^ < 10.4. For A^ > 10.4, they exhibit collapse, again through the splitting into a 
set of two S = solitons which then blow up intrinsically. Lastly, we have also checked that 
the zero-vorticity sohtons in Eq. (jH)) exhibit collapse exactly when A^ > A^i^x^'* = 5.85, as 
was assumed above. 

At a late stage of the collapse, when the condensate will shrink to a size ~ a^, the above 
derivation of the 2D equation (jSJ from its 3D counterpart Q will break down, and the 
initially quasi-2D collapse will switch into a faster 3D blow-up mode. Direct simulations of 
the full 3D equation pi|l corroborate this expectation, as will be reported in detail elsewhere. 



12 



V. VORTEX SOLITONS IN THE THREE-DIMENSIONAL GROSS-PITAEVSKII 
EQUATION 



In this section, we aim to compare the above results, obtained in the framework of the 2D 
approximation, with what can be found from the full 3D equation where the quasi-2D 
case corresponds to large Q. Systematic presentation of the results for the 3D model will be 
a subject of a separate work, while here we report, in a brief form, findings which are most 
relevant for the comparison with the 2D approximation. 

Stationary solutions for the 3D solitons with embedded vorticity S and chemical potential 
fi are looked for in the form [cf. Eq. (jl4j) ] \1/ = R{r, z) exp {iS6 — ifit), where r,z,d are 
cylindrical coordinates, and the real function R obeys the equation 

^ + '-f + ^+(2,-^-r^-n^AR-2R^ = 0, (26) 
ar^ r or oz'^ \ / 

supplemented by obvious boundary conditions: R{r,z) when r ^ oo or |2;| ^ oo, 

and i? ~ r"^ for r ^ at fixed z (we again define S* to be > 0). After finding a family 

of the vortices from numerical solution of Eq. ()26|). their stability was investigated via the 

computation of eigenvalues for small perturbations. 

First of all, the linear version of Eq. can be considered similar to how it was done 

above for the 2D case: the solutions are built as linear combinations of products of the wave 

functions of the harmonic oscillator, cf. Eq. (fTBj) : 

^,ki{x,y,z,t) = e-'^^'+'+^^+'/'^''H,{x)My)M^z), (27) 

where / = 0, 1, 2, ... is the quantum number in the 2;-direction. Being interested here in the 
"flattest" states (ones closest to the pancake configuration in the 3D space), we will set 
/ = 0. Then, wave functions with given vorticity S in the {x,y) plane are constructed as 
combinations of factorized wave functions (j27jl . subject to the same constraint as in the 2D 
case, i.e., j + k = S . From this restriction and from the form of the chemical potential in 
the linear-limit solution (j27|) with 1 = 0, fi = j + k + l + Q/2,it follows that the chemical 
potential of solutions to the full nonlinear equation (pUj) . obtained as a continuation starting 
from the linear limit, obeys a limitation 

1^ < /^max = S + 1 + n/2, (28) 

cf. Eq. (1201). 
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Numerically found families of the 3D soliton solutions with S = 1 are characterized by 
yu(A^) and E{N) curves displayed in Fig. [71 one for the anisotropic model, with Q = 10, 
and one for its isotropic counterpart, with Q = 1. In the former case [which, as a matter of 
fact, is not a strongly anisotropic one, as the corresponding aspect ratio of the confinement 
lengths in the (x, y) and z directions is vTO, pursuant to Eq. ((T))], the existence and stability 
intervals of the fundamental vortices have been found to be, respectively, < < 14.46, 
and 

0<iV<iV(3^)|f,=io = 5.95 (29) 

[cf. Eqs. fl23j) and (j^^ in the 2D model], hence the relative size of the stability area is 
^ 0.41. This should be compared to the relative size of the stability area in the 2D limit 
reported above, which is ~ 1/3. For the isotropic model {Q = 1), the existence and stability 
intervals are, respectively, < < 23, and 

0<iV< Ar|^)|o=i^l5, (30) 

their ratio being ^ 0.65. 

Thus, we conclude that the 2D limit yields a reasonable approximation for the pancake- 
shaped configurations (the relative difference between its prediction and the direct numerical 
result for the moderately anisotropic model with the aspect ratio a/TO is ~ 20%). Another 
noteworthy conclusion is that the relative size of the stability area essentially increases as 
the 3D model approaches the isotropic configuration; in the isotropic limit, the relative size 
of the stability region is almost twice as large as in the 2D limit. In terms of the maximum 
value of the norm, A^max? up to which the vortex soliton remains stable, the difference is still 
bigger: A^max for Q = 1 exceeds its counterpart for = 10 by a factor of ~ 2.5. 

To further illustrate the effect of the degree of asymmetry on the stability of vortex 
solitons, in Fig. |Hlwe display the spectrum of unstable eigenvalues for the 3D solitons with 
= 1, in the same two CcLSGS cLS considered above, i.e., the moderately anisotropic one with 
Q = 10, and isotropic with = 1. Comparison of this figure with the respective picture 
from Fig. |21^a) for the 2D model shows a gradual deformation of the spectrum with the 
transition from the 2D limit to the isotropic 3D configuration. An important observation is 
that the transition from the 2D approximation to the full 3D model does not lead to any new 
instability eigenmode (at least, up to the isotropic limit, Q = 1); incidentally, the bending 
instability of the vortex core does not occur either. 
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As well as in the 2D case, vortex solitons with S* > 2 [at least, the "flattest" ones, 
corresponding to Z = in the linear limit, see Eq. (jTTj) ] are all unstable. As concerns the 
fundamental solitons with S* = 0, we could easily reproduce a known result stating that, as 
well as in the 2D model, they are stable within their full existence region, which is < 3.8 
for Q = 10, and A^ < 7.3 for Q = 1. Comparing these values with ones in Eqs. (|29|) and 
fl3U|) ■ we conclude that the introduction of the vorticity enhances the stability limit by a 
factor of 1.6 for Q = 10, and by 2 for = 1. Recall that, in the 2D model, the same 
enhancement factor was 1.3 [if its definition does not include the extra region ()25|) of the 
quasi-stable dynamical states]; this again demonstrates a general trend to the expansion of 
stability regions with the transition from the 2D limit to the fully isotropic 3D case. 



VI. CONCLUSION 



We have reported results of accurate investigation of the stability of localized vortices in 
Bose-Einstein condensates with self-attraction, trapped in the nearly two-dimensional con- 
figuration. We have briefiy recapitulated the derivation of the effective 2D Gross-Pitaevskii 
equation from its full 3D counterpart, and demonstrated that, in the 2D limit, the vortices 
with S = 1 are stable against azimuthal splitting in a third of their existence region, i.e., for 
< A^ < (l/3)A^i^x^'*, in terms of the solution's norm (which is proportional to the number 
of atoms trapped in the vortex soliton). In the adjacent interval, (l/3)A'max < N < 
0.43A'max , the unstable vortex does not collapse, but rather demonstrates stable quasi- 
periodic dynamics, featuring cycles of splitting into two fragments and recombination, which 

(S=l) 

is a novel dynamical regime. For A^ > 0.43A^max , the two segments produced by the split- 
ting fail to recombine, as they collapse intrinsically. The stable vortex soliton created in the 
pancake-shaped condensate of ^Li may contain up to 10, 000 atoms. 

The results were checked against direct simulations of the 2D Gross-Pitaevskii equation, 
and, which is crucially important for the verification of their physical relevance, they were 
compared with the stability analysis of localized vortices in the full 3D equation. It was 
concluded that, for the moderately anisotropic 3D configuration with the aspect ratio v^TO ~ 
3.2, the stability interval of the 5 = 1 vortices occupies ~ 40% of their existence region 
(hence, the 2D limit provides for quite a reasonable approximation). In the 3D isotropic 
limit, the stability interval expands to 65% of the existence domain. The enhancement factor 
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for the collapse threshold, caused by the introduction of vorticity S — 1, was found too; in 
the 3D isotropic configuration it attains the value of 2. All higher-order localized vortices, 
with S > 2, are completely unstable, in the 2D limit and in the full 3D Gross-Pitaevskii 
equation alike. 
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FIG. 1: Chemical potential /x (a) and energy E (b) versus normalized number of atoms N for 
solutions with vorticities S = 0, 1, and 2, as found from Eq. H15|) with g = —1. Solid and dashed 
lines show stable and unstable branches of the solutions (see below). The arrow indicates the point 
where the 5 = 1 vortices loose their stability. 

FIG. 2: (a) The real part of the perturbation growth rate A for different azimuthal indices L. (b) 
The real and imaginary parts (solid and dashed lines, respectively) of the eigenvalues responsible 
for the Hamiltonian-Hopf bifurcation which destabilizes the vortices with S" = 1 at the critical 
point, /i = /icr ~ 1.276. This value and the bifurcation point are marked by arrows. 

FIG. 3: Grey-scale plots illustrating recovery of the perturbed stable vortex with S = 1 for fi = 1.4. 
(a,b): Intensity and phase distributions in the initial configuration including random noise. (c,d): 
The same in the self-cleaned vortex soliton at t = 120. The norms of the input and output solitons 
are = 6.641 and = 6.629, respectively (the small loss of the norm is due to absorbers placed 
at borders of the integration domain). 

FIG. 4: Formation of a stable vortex with S = 1 from a Gaussian input. (a,b): The intensity 
and phase distributions in the initial Gaussian with a nested vortex. (c,d): The intensity and 
phase distributions in the output {t = 700). The input and output configurations have the norm, 
respectively, N = 6.406 and N = 6.392. 

FIG. 5: Quasi-periodic evolution of the amplitude of the unstable S = 1 vortex for fi = 1.2. The 
labels a, b, c and d point at values of the amplitude at t = 0, 100, 140, and 180 (see Fig. 6). 

FIG. 6: Regular evolution of the unstable vortex with S = 1 for fi = 1.2 initiated by random noise, 
which shows periodic splitting into two fragments and their subsequent recombination: (a) t = 0, 
(b) t = 100, (c) t = 140, and (d) t = 180. Here, values of the input and output norm are = 8.476 
(at t = 0) and iV = 8.460 (at t = 240). 

FIG. 7: The chemical potential (a) and energy (b) for families of three-dimensional solitons with 
the embedded vorticity S" = 1 vs. the normalized number of atoms, with the moderately anisotropic 
($7 = 10) and isotropic {0, = 1) confining potentials. The continuous and dashed parts of the curves 
show, respectively, stable and unstable parts of the families. 
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FIG. 8: The same as in Fig. Efa), but for vortex solitons with 5" = 1 in the full three-dimensional 
model: (a) O = 10; (b) 0, = 1. The arrows show the edge of the existence regions of the solitons 
on the scale of the chemical potential, as per Eq. (|28j) . 
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